library(broom)
library(dplyr)
library(purrr)
library(tidyr)
detach("package:plyr", unload = TRUE)

df3_het <- df3_het %>%
  mutate(profileBI = factor(profileBI))

controls <- c("gender", "age", "degree", "ethnicity", "partyID", "partners", "rightwing", "attentionpass")
control_labels <- c(
gender = "Gender",
age = "Age",
degree = "Education",
ethnicity = "Ethnicity",
partyID = "Partisan ID",
partners = "N. of partners",
rightwing = "Ideology",
attentionpass = "Attention checks")


control_sets <- unlist(lapply(1:length(controls), function(i) {
  combn(controls, i, simplify = FALSE)
}), recursive = FALSE)

spec_results <- tibble(control_set = seq_along(control_sets)) %>%
  mutate(
    controls = map(control_set, ~ control_sets[[.x]]),
    formula_chr = map_chr(controls, ~ paste("profile_recommend ~", paste(c("profileBI", .x), collapse = " + "))),
    formula_obj = map(formula_chr, as.formula),
    model = map(formula_obj, ~ lm(.x, data = df3_het)),
    tidy = map(model, ~ broom::tidy(.x))
  ) %>%
  unnest(tidy) %>%
  filter(term == "profileBI2") %>%
  mutate(
    lower90 = estimate - 1.645 * std.error,
    upper90 = estimate + 1.645 * std.error,
    lower95 = estimate - 1.96 * std.error,
    upper95 = estimate + 1.96 * std.error,
    sig_05 = p.value < 0.05,
    model_id = row_number()
  ) %>%
  arrange(estimate) %>%
  mutate(model_id = row_number())

spec_results <- spec_results %>%
  arrange(estimate) %>%
  mutate(model_id = row_number())

min_estimate <- min(spec_results$estimate, na.rm = TRUE)
median_estimate <- median(spec_results$estimate, na.rm = TRUE)
max_estimate <- max(spec_results$estimate, na.rm = TRUE)


spec_matrix <- spec_results %>%
  select(model_id, controls) %>%
  mutate(
    control_flags = map(controls, ~ set_names(rep(1, length(.x)), .x))  
  ) %>%
  select(model_id, control_flags) %>%
  unnest_wider(control_flags, names_repair = "unique") %>%  
  pivot_longer(
    cols = -model_id,
    names_to = "control",
    values_to = "included"
  ) %>%
  mutate(
    included = factor(ifelse(is.na(included), 0, included))  )
min_estimate <- min(spec_results$estimate, na.rm = TRUE)
median_estimate <- median(spec_results$estimate, na.rm = TRUE)
max_estimate <- max(spec_results$estimate, na.rm = TRUE)


spec_matrix <- spec_results %>%
  select(model_id, controls) %>%
  mutate(
    control_flags = map(controls, ~ set_names(rep(1, length(.x)), .x))  
  ) %>%
  select(model_id, control_flags) %>%
  unnest_wider(control_flags, names_repair = "unique") %>%  
  pivot_longer(
    cols = -model_id,
    names_to = "control",
    values_to = "included"
  ) %>%
  mutate(
    included = factor(ifelse(is.na(included), 0, included))  )

p_curve <- ggplot(spec_results, aes(x = model_id, y = estimate)) +
  geom_hline(yintercept = median_estimate, linetype = "solid", linewidth=1, color = "#660033") +
  geom_linerange(aes(ymin = lower95, ymax = upper95), alpha = 0.5, linewidth = 1, color = "#E0115F") +
  geom_linerange(aes(ymin = lower90, ymax = upper90), alpha = 0.9, linewidth = 1, color = "#E0115F") +
  geom_point(size = 2, fill="white", color = "#E0115F", pch=21) +  
  geom_hline(yintercept = 0, linetype = "dashed", color = "#660033") +
  labs(subtitle="Effect on 'would encourage a friend to date'",
    title = "255 different specifications with a median ATE of -8.9 percentage-points",
    x = "",
    y = "Average treatment effect",
    shape = "p < 0.05") +
  annotate(geom="label", x = 50, y =(median_estimate+.01), size = 3, color = "#660033", fontface=2,
           label = "Median point-estimate: -0.089")+
  theme_ggdist()+
  theme(
    plot.title = element_text(size=14, face="bold"))

spec_matrix <- spec_matrix %>%
  mutate(control = recode(control, !!!control_labels))

p_controls <- ggplot(spec_matrix, aes(x = model_id, y = fct_rev(factor(control)), fill = included)) +
  geom_tile(color = "white", linewidth = 0.2) +
  scale_fill_manual(values = c("0" = "gray90", "1" = "#E0115F")) +
  labs(x = NULL, y = "Included Controls") +
  theme_ggdist()+
  theme(
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    panel.grid = element_blank(),
    plot.title = element_text(size=14, face="bold"),
    legend.position = "none")

p_curve / p_controls +
  plot_layout(heights = c(3, 1))+
  plot_annotation(title = 'Study 3: Multiverse analysis I',
                  caption = "Confidence intervals at 95% & 90%",
                  theme = theme(plot.title = element_text(face = "bold", size=18),
                                plot.subtitle = element_text(size=14)))
ggsave("figures/appendix/FigureC6.png", width=8.5)


spec_results <- tibble(control_set = seq_along(control_sets)) %>%
  mutate(
    controls = map(control_set, ~ control_sets[[.x]]),
    formula_chr = map_chr(controls, ~ paste("profile_woulddate ~", paste(c("profileBI", .x), collapse = " + "))),
    formula_obj = map(formula_chr, as.formula),
    model = map(formula_obj, ~ lm(.x, data = df3_het)),
    tidy = map(model, ~ broom::tidy(.x))
  ) %>%
  unnest(tidy) %>%
  filter(term == "profileBI2") %>%
  mutate(
    lower90 = estimate - 1.645 * std.error,
    upper90 = estimate + 1.645 * std.error,
    lower95 = estimate - 1.96 * std.error,
    upper95 = estimate + 1.96 * std.error,
    sig_05 = p.value < 0.05,
    model_id = row_number()
  ) %>%
  arrange(estimate) %>%
  mutate(model_id = row_number())

spec_results <- spec_results %>%
  arrange(estimate) %>%
  mutate(model_id = row_number())

min_estimate <- min(spec_results$estimate, na.rm = TRUE)
median_estimate <- median(spec_results$estimate, na.rm = TRUE)
max_estimate <- max(spec_results$estimate, na.rm = TRUE)


spec_matrix <- spec_results %>%
  select(model_id, controls) %>%
  mutate(
    control_flags = map(controls, ~ set_names(rep(1, length(.x)), .x))  
  ) %>%
  select(model_id, control_flags) %>%
  unnest_wider(control_flags, names_repair = "unique") %>%  
  pivot_longer(
    cols = -model_id,
    names_to = "control",
    values_to = "included"
  ) %>%
  mutate(
    included = factor(ifelse(is.na(included), 0, included))  )
min_estimate <- min(spec_results$estimate, na.rm = TRUE)
median_estimate <- median(spec_results$estimate, na.rm = TRUE)
max_estimate <- max(spec_results$estimate, na.rm = TRUE)


spec_matrix <- spec_results %>%
  select(model_id, controls) %>%
  mutate(
    control_flags = map(controls, ~ set_names(rep(1, length(.x)), .x))  
  ) %>%
  select(model_id, control_flags) %>%
  unnest_wider(control_flags, names_repair = "unique") %>%  
  pivot_longer(
    cols = -model_id,
    names_to = "control",
    values_to = "included"
  ) %>%
  mutate(
    included = factor(ifelse(is.na(included), 0, included))  )

p_curve <- ggplot(spec_results, aes(x = model_id, y = estimate)) +
  geom_hline(yintercept = median_estimate, linetype = "solid", linewidth=1, color = "#660033") +
  geom_linerange(aes(ymin = lower95, ymax = upper95), alpha = 0.5, linewidth = 1, color = "#E0115F") +
  geom_linerange(aes(ymin = lower90, ymax = upper90), alpha = 0.9, linewidth = 1, color = "#E0115F") +
  geom_point(size = 2, fill="white", color = "#E0115F", pch=21) +  
  geom_hline(yintercept = 0, linetype = "dashed", color = "#660033") +
  labs(subtitle="Effect on 'would personally date'",
       title = "255 different specifications with a median ATE of -11 percentage-points",
       x = "",
       y = "Average treatment effect",
       shape = "p < 0.05") +
  annotate(geom="label", x = 50, y =(median_estimate+.01), size = 3, color = "#660033", fontface=2,
           label = "Median point-estimate: -0.105")+
  theme_ggdist()+
  theme(
    plot.title = element_text(size=14, face="bold"))

spec_matrix <- spec_matrix %>%
  mutate(control = recode(control, !!!control_labels))

p_controls <- ggplot(spec_matrix, aes(x = model_id, y = fct_rev(factor(control)), fill = included)) +
  geom_tile(color = "white", linewidth = 0.2) +
  scale_fill_manual(values = c("0" = "gray90", "1" = "#E0115F")) +
  labs(x = NULL, y = "Included Controls") +
  theme_ggdist()+
  theme(
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    panel.grid = element_blank(),
    plot.title = element_text(size=14, face="bold"),
    legend.position = "none")

p_curve / p_controls +
  plot_layout(heights = c(3, 1))+
  plot_annotation(title = 'Study 3: Multiverse analysis II',
                  caption = "Confidence intervals at 95% & 90%",
                  theme = theme(plot.title = element_text(face = "bold", size=18),
                                plot.subtitle = element_text(size=14)))
ggsave("figures/appendix/FigureC7.png", width=8.5)

